************************************************************
****************** REPLICATION CODE ************************
************************************************************
* For: Rigterink, The Wane of Command, APSR				   *
* Simulations of expanded sample size					   *
************************************************************
* Content of file:										   *
/*														   *
1. Preliminaries and setting directory					   *
2. Finding distribution of terrorist attacks			   *
3. Finding distribution of drone strikes				   *
4. Main results											   *
5. Setting parameters for each terrorist group 			   *
6. Simulation of expanded sample size					   *
7. Graph of results										  */
************************************************************

/*----------------------------------------------------------
---------- Preliminaries and setting directory -------------
----------------------------------------------------------*/

clear all
version 14
set more off

*******************************
*** Set your directory here ***
*******************************
cd "[directory]"

cap log close
log using log\Rigterink_drones_replication_appendix_samplesize.log, replace

use "dta/Rigterink_drones_replication_maindataset.dta", clear


/*----------------------------------------------------------
------- Finding distribution of terrorist attacks ----------
----------------------------------------------------------*/

* Test Poisson against NBREG and ZINB

forvalues i = 1(1)13 {
	
	cap poisson terratt_unlog if groupid==`i', iterate(250)
	scalar ll_pois`i' = e(ll)

	cap nbreg terratt_unlog if groupid==`i', iterate(250)
	scalar ll_nb`i' = e(ll)
		
	cap zinb terratt_unlob if groupid==`i', inflate(_cons) , iterate(250)
	scalar ll_zinb`i' = e(ll)
	
	di "NB vs ZINB"
	di "LL NB `i': " ll_nb`i'
	di "LL ZINB `i': " ll_zinb`i'
	di "Prob > chi2 = " chi2tail(2, 2*(abs(ll_zinb`i'-ll_nb`i')))
	
	di "NB vs POISSON"
	di "LL NB `i': " ll_nb`i'
	di "LL POISSON `i': " ll_pois`i'
	di "Prob > chi2 = " chi2tail(2, 2*(abs(ll_pois`i'-ll_nb`i')))
	
	}

*** No evidence that the ZINB performs any better than the NB
* In fact, NB appears to perform slightly better
*** Evidence that NB performs better than Poisson


/*----------------------------------------------------------
--------- Finding distribution of drone strikes ------------
----------------------------------------------------------*/

* Test Poisson against NBREG and ZINB

forvalues i = 1(1)13 {
	
	cap poisson dronestrike if groupid==`i', iterate(250)
	scalar ll_pois`i' = e(ll)

	cap nbreg dronestrike if groupid==`i', iterate(250)
	scalar ll_nb`i' = e(ll)
		
	cap zinb dronestrike if groupid==`i', inflate(_cons) , iterate(250)
	scalar ll_zinb`i' = e(ll)
	
	di "NB vs ZINB"
	di "LL NB `i': " ll_nb`i'
	di "LL ZINB `i': " ll_zinb`i'
	di "Prob > chi2 = " chi2tail(2, 2*(abs(ll_zinb`i'-ll_nb`i')))
	
	di "NB vs POISSON"
	di "LL NB `i': " ll_nb`i'
	di "LL POISSON `i': " ll_pois`i'
	di "Prob > chi2 = " chi2tail(2, 2*(abs(ll_pois`i'-ll_nb`i')))
	
	}

*** No evidence that the ZINB performs any better than the NB
*** Evidence that NB performs better than Poisson for 5/13 groups


/*----------------------------------------------------------
---------------------- Main results ------------------------
----------------------------------------------------------*/

* Generating results in main paper

reghdfe terratt L(0/6).bigfishtarget L(0/6).bigfishdied L(0/6).dronestrike F(1/6).bigfishtarget F(1/6).bigfishdied F(1/6).dronestrike, absorb(i.groupid i.periodid) vce(, bw(12))
matrix realb = e(b)
matrix list realb


/*----------------------------------------------------------
------ Setting parameters for each terrorist group ---------
----------------------------------------------------------*/

* Matrix for parameters
matrix para = J(13, 11, .)

* Variable for miss
gen miss = 1 if bigfishtarget==1 & bigfishdied==0
replace miss = 0 if miss!=1

* Setting parameters and populating matrix
forvalues i=1(1)13 {

	qui nbreg terratt_unlog if groupid==`i', iterate(100)
	matrix para[`i', 1]= e(b)
	matrix para[`i', 3]= e(alpha)
	matrix para[`i', 4]= 1/para[`i', 3]
	matrix para[`i', 5]= exp(para[`i', 1])
	
	qui sum bigfishdied if groupid==`i'
	matrix para[`i', 6]= r(mean)
	qui sum miss if groupid==`i'
	matrix para[`i', 7]= r(mean)
	
	qui nbreg dronestrike if groupid==`i', iterate(100)
	matrix para[`i', 8]= e(b)
	matrix para[`i', 9]= e(alpha)
	matrix para[`i', 10]= 1/para[`i', 9]
	matrix para[`i', 11]= exp(para[`i', 8])
	
}


/*----------------------------------------------------------
----------- Simulation of expanded sample size -------------
----------------------------------------------------------*/

set matsize 1000

foreach Y in 10 25 50 {
	
	local count = `count'+1
	
	di "Additional years: " `Y'
	
	local obs = 1733+`Y'*12*13

	set obs `obs'

	qui replace groupid = (_n-_N) - round((_n-_N)-7, 13) if groupid==.
	qui replace periodid = floor((_n+139+12)/13) if periodid==.

	forvalues R = 1(1)1000 {
		
		di "Round: " `R'
		
		local seed = 19860223+`R'
		set seed `seed'
	
		qui gen simdep = .
		qui gen simdied = .
		qui gen simtarget = .
		qui gen simdronestrike = . 
	
		forvalues i = 1(1)13 {
			qui gen xg = rgamma(para[`i', 4], para[`i', 3]) if groupid==`i'
			qui gen xbg = xg*para[`i', 5] if groupid==`i'
			qui replace simdep = rpoisson(xbg) if groupid==`i'
	
			qui drop xg xbg
	
			qui replace simdied = rbinomial(1, para[`i', 6]) if groupid==`i'
			qui replace simtarget = 1 if simdied==1
			qui replace simtarget = rbinomial(1, para[`i', 7]) if simtarget==. & groupid==`i'
		
			qui gen xg = rgamma(para[`i', 10], para[`i', 9]) if groupid==`i'
			qui gen xbg = xg*para[`i', 11] if groupid==`i'
			qui replace simdronestrike = rpoisson(xbg) if groupid==`i'
			
			qui drop xg xbg
		
		}

		qui replace simdep = 0 if simdep==.
		qui replace simdied = 0 if simdied==.
		qui replace simtarget = 0 if simtarget==.
		qui replace simdronestrike = 0 if simdronestrike==.
		
		qui replace simdep = ln(simdep+1)
		
		qui xtset groupid periodid

		qui reghdfe simdep L(0/6).simtarget L(0/6).simdied L(0/6).simdronestrike F(1/6).simtarget F(1/6).simdied F(1/6).simdronestrike, absorb(i.groupid i.periodid) vce(, bw(12))
		matrix b=e(b)
	
		if `R'==1 {
			matrix simb`Y' = b
		}
		else {
			matrix simb`Y' = (simb`Y' \ b)
		}
	
		qui drop sim*
	}
	
}

* Storing results from matrices in separate dataset

clear

foreach Y in 10 25 50 {

	svmat simb`Y'
	rename simb`Y'* simb*
	gen ssize = `Y'
	
	if `Y'==10 {
		save dta/Rigterink_drones_replication_samplesize.dta, replace
		clear
	}
	
	else {
		append using dta/Rigterink_drones_replication_samplesize.dta
		save dta/Rigterink_drones_replication_samplesize.dta, replace
		clear
	}
}


clear
svmat realb

rename realb* simb*
gen ssize = 99 

append using dta/Rigterink_drones_replication_samplesize.dta
save dta/Rigterink_drones_replication_samplesize.dta, replace


/*----------------------------------------------------------
------------------- Graph of results -----------------------
----------------------------------------------------------*/

use dta/Rigterink_drones_replication_samplesize.dta, clear

for var simb*: gen p5_X = X \ rename X p95_X

collapse (p5) p5* (p95) p95*, by(ssize)

reshape long p5_simb p95_simb, i(ssize) j(time)

keep if (time>27 & time<34) | (time>7 & time<15)

replace time = -1*(time-27) if time>20
replace time = time-8 if time>7

sort ssize time

label define time -6 "t-6" -5 "t-5" -4 "t-4" -3 "t-3" -2 "t-2" -1 "t-1" 0 "t" 1 "t+1" 2 "t+2" 3 "t+3" 4 "t+4" 5 "t+5" 6 "t+6"
label values time time

graph twoway ///
(line p5_simb time if ssize==10, lcolor(black) lpattern(solid)) || ///
(line p95_simb time if ssize==10, lcolor(black) lpattern(solid)) || ///
(line p5_simb time if ssize==25, lcolor(black) lpattern(dot)) || ///
(line p95_simb time if ssize==25, lcolor(black) lpattern(dot)) || ///
(line p5_simb time if ssize==50, lcolor(black) lpattern(dash)) || ///
(line p95_simb time if ssize==50, lcolor(black) lpattern(dash)) || ///
(scatter p95_simb time if ssize==99, mcolor(black) msymbol(D)), ///
plotregion(style(none)) graphregion(ifcolor(white) fcolor(white) color(white) icolor(white)) ///
legend(order(1 "90% CI, 10 additional years" 3 "90% CI, 25 additional years" 5 "90% CI, 50 additional years" 7 "Actual coefficient estimates")) ///
title("Sample size simulations") ytitle("Coefficient on hit") xtitle("") ///
xlabel(-6(1)6, valuelabel)

graph export graphs/Rigterink_AFigB8.pdf, replace


*******************
*** END OF FILE ***
*******************
